arXiv: 1507.06764v2 [physics.comp-ph] 4 Feb 2016 


Sampling Free Energy Surfaces as Slices by Combining 
Umbrella Sampling and Metadynamics 

Shalini Awasthi, Venkat Kapil and Nisanth N. Nair^ 

Department of Chemistry 

Indian Institute of Technology Kanpur, 208016 Kanpur, India 

Keywords: Metadynamics, Umbrella Sampling, Reweighting, Weighted Histogram 
Analysis, Free energy calculations 


1* 


Corresponding Author: nnair@iitk.ac.in 



Abstract 


Metadynamics (MTD) is a very powerful technique to sample high-dimensional free 
energy landscapes, and due to its self-guiding property, the method has been success¬ 
ful in studying complex reactions and conformational changes. MTD sampling is based 
on filling the free energy basins by biasing potentials and thus for cases with flat, broad 
and unbound free energy wells, the computational time to sample them becomes very 
large. To alleviate this problem, we combine the standard Umbrella Sampling (US) 
technique with MTD to sample orthogonal collective variables (CVs) in a simultane¬ 
ous way. Within this scheme, we construct the equilibrium distribution of CVs from 
biased distributions obtained from independent MTD simulations with umbrella po¬ 
tentials. Reweighting is carried out by a procedure that combines US reweighting and 
Tiwary-Parrinello MTD reweighting within the Weighted Histogram Analysis Method 
(WHAM). The approach is ideal for a controlled sampling of a CV in a MTD simulation, 
making it computationally efficient in sampling flat, broad and unbound free energy 
surfaces. This technique also allows for a distributed sampling of a high-dimensional 
free energy surface, further increasing the computational efficiency in sampling. We 
demonstrate the application of this technique in sampling high-dimensional surface 
for various chemical reactions using ab initio and QM/MM hybrid molecular dynam¬ 
ics simulations. Further, in order to carry out MTD bias reweighting for computing 
forward reaction barriers in ab initio or QM/MM simulations, we propose a computa¬ 
tionally affordable approach that does not require recrossing trajectories. 
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1 Introduction 


Metadynamics (MTD) is a very powerful tool to sample complex free energy landscapes 
of complex chemical reactions, phase transitions and conformational changes. 0 § 
MTD is a biased sampling approach where sampling of a selected set of collective 
variables (CVs) is enhanced by introducing slowly grown smooth history dependent 
repulsive potentials along the trajectory of CVs. In MTD, the underlying potential 
energy of a system li(R) is modified to li(R) + V b (s, t), where V b (s, t) is the biasing 
potential applied along the CVs, s = s(R), at any instance of the simulation t. Typically, 
V b (s, t) is constructed by the sum of spherical Gaussians deposited over time, as 


V b (s,f) = zv T exp 


T<t 


2(5s) 2 


(1) 


Here w T and 5s are the height and the width parameters defining the Gaussian potential 
added at some time r. In Well-Tempered MTD (WT-MTD) 

V h (s,t ) 


w T = coqTq exp 


/c b AT 


( 2 ) 


where coq is the initial rate of deposition of the bias, tq is the time step at which Gaussian 
potentials are augmented, and AT is a parameter. The advantage of WT-MTD is that a 
systematic convergence in free energy can be achieved - at the limit t —> oo, the biasing 
potential V b (s, f) does not vary much, and 

AT 


limy (8,0+/'— f(s) 


(3) 


where f is a constant. Thus, a converged free energy surface can be constructed as 


T(s) = —tx lim V b (s, t) + f (4) 

f—S>oo 

where a = (T + AT) / AT and / is some other constant. [7] 

In MTD simulations the total simulation time required to sample a free energy land¬ 
scape depends exponentially on the number of CVs. Larger the volume of the free 
energy wells, more is the time required to fill them and to see transitions from one well 
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to the other. Generally, MTD simulations are carried out with 2 or 3 CVs [3]. Although, 
most of the reactions or structural changes can be sampled using 2 or 3 CVs, one finds 
several other orthogonal coordinates which have hidden barriers, leading to serious er¬ 
rors in the free energies and poor convergence. Thus inclusion of a few more CVs in 
sampling would help to accelerate sampling of orthogonal coordinates. In this spirit, 
MTD has been combined with parallel tempering method |8J, and a technique called 
bias-exchange MTD |j9 10 [ has been introduced. 

An ideal technique to compute free energy along a known reaction path is the Um¬ 
brella Sampling (US) method, |TT| where a number of time independent harmonic 
basing potentials are applied to obtain biased probability distribution of CVs. Here the 
umbrella biasing potential placed at s/„ 


Wj(s) = ^ K h(s~Sh) Z , h — 1,- ■ ■ ,M 


(5) 


is used to construct a set of biased probability distributions [Pj ,(s)}, from M inde¬ 
pendent MD simulations with varying biases. Subsequently, {P/ 7 (s)} is reweighted to 
obtain the equilibrium probability distribution {Ph( s)} for all the M windows, and are 
subsequently combined to get the total distribution P(s) through the Weighted His¬ 


togram Analysis Method (WHAM). [12,13] The equilibrium free energy surface is then 
constructed by 


P(s) = -^lnP(s)+/ (6) 

where /3 = and / is some arbitrary constant. 

In Figure [lj, we sketch some practical limitations of standard MTD simulation in 
sampling flat potentials. For potentials like those shown in Figure [l^ and b, MTD 
simulations fail to sample transition from reactant to product wells. For the case in 
Figure |lj:, a large number of Gaussian potentials has to be filled in the reactant basin 
due to its broadness. Chemical reactions in systems such as weakly bound Michaelis 
complexes in enzymes, weakly bound molecular complexes, and in general A+B type 
reactions in solutions and gasphase have such free energy topology. Here, MTD spends 
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most of the computational time in filling uninteresting parts of the free energy wells. 
Partly, these problems can be circumvented by using repulsive wall potentials. How¬ 
ever, wall potentials introduce boundary effects and require corrections for the bound¬ 
ary effects. fll4|[T5fl Often it is difficult to scrutinize whether the minimum observed 
near the boundary is due to the artifact of the wall potential or actual. Also, the co¬ 
ordinate that requires controlled sampling using wall potentials has to be defined as 
CVs (if they are not part of the CVs already). Alternatively, the coordinate along which 
the free energy is flat or broad can be efficiently sampled using US since the sampling 
range of the coordinate can be controlled in US. However, to sample bound orthogonal 
coordinates with hidden barriers, MTD is ideal because of its self-guiding nature. Thus 
for an efficient sampling of a high dimensional free energy landscape that has features 
as in Figure [I] along certain CVs, a combination of US with MTD, where these methods 
simultaneously sample orthogonal coordinates, is ideal. 

Here we report a procedure to carry out such hybrid simulations and to obtain 
free energy landscape in the full CV space. We name this technique of combining US 
and MTD as Well-Sliced MTD (WS-MTD). Here M independent MTD simulations are 
carried out with the bias 
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to sample the CV space 



s = 


Here W^(s a ) and V^(sp, t) are given by Equation ([5]) and Equation (|T|), respectively. M 
number of umbrella biases are placed along the s a coordinates, and for each of these 
umbrella h, we carry out MTD simulation sampling the coordinates. Each US+MTD 
simulation samples a slice of the high-dimensional free energy surface. Subsequently 
we combine the biased probability distributions from M different US+MTD simulations 
to construct the total equilibrium distribution and the free energy surface. 

Reweighting is not straightforward for MTD, as the biasing potential is time depen¬ 
dent and the sampling weights change with simulation time. Tiana [16] proposed a 
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strategy to obtain ensemble averages from MTD based on a time-dependent reweight¬ 
ing scheme. Employing the convergence property of WT-MTD, a more systematic strat¬ 


egy was reported by Bonomi et al. [ 171 Laio and co-workers have put forward a differ¬ 
ent reweighting scheme in the framework of bias-exchange MTD. [ |T0| Recently, a very 


simple and efficient reweighting scheme was reported by Tiwary and Parrinello [18] 
based on time dependent weights directly computable from WT-MTD simulations. 

Combining US with MTD has been reported first by Frilizola and co-workers 19 [ 
where a MTD reweighting proposed by Bonomi et al. |17| was used. US technique 
has been combined with various sampling techniques; see Ref. [201 for a review and 


Ref. 121 [ for a recent example. There were also attempts to do US corrections on the 
free energy surface obtained from MTD. [ 


25[ 


Here we combine US and MTD to sample orthogonal coordinates and reconstruct 
the free energy surface by a combination of Tiwary-Parrinello reweighting scheme and 
US reweighting within WHAM. First we carefully study the method on a model poten¬ 
tial for which the exact free energy barriers are known. Then we extend our study to 
various problems that are part of the ongoing research in our laboratory, where normal 
MTD has failed or has poor performance in sampling high dimensional free energy 
landscape due to the problems shown in Figure [l| In this respect, we first perform ab 
initio MD simulation of formation of cyclobutene from 1,3-butadiene using WS-MTD 
and compare its performance with WT-MTD. This serves as an ideal example as its 
free energy landscape has very broad and deep reactant basin. We then applied WS- 
MTD to model ligand exchange reactions of Pd complex in aqueous solution where the 
standard MTD simulations are known to fail. [261 This example shows how controlled 
sampling of a high dimensional free energy landscape can be achieved by WS-MTD. 
Finally, we demonstrate the efficiency of WS-MTD in modeling an enzymatic reaction 
using QM/MM method. Coordination of water molecule to one of the catalytic Zn ions 
in the active site of New Delhi metallo jS-lactamase 1 (NDM-1) |27j is modeled, for 
which our earlier attempts using WT-MTD and non-tempered MTD simulations were 
unsuccessful. 
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2 Methods and Models 


2.1 Reweighting Scheme in WS-MTD 

Here we discuss how to obtain the reweighted distribution from WS-MTD simulations. 
Let us consider a problem where we are interested in computing the free energy surface 
F(si,S 2 ) by sampling the CVs Si and s 2 . Consider that, the coordinate Si is sampled 
using US, while S 2 is simultaneously sampled by a one-dimensional WT-MTD. If M 
umbrella potentials are placed along the si coordinate, we carry out M independent 
MTD simulations, each using the Hamiltonian 

H / 7 s (R,P) = H°(R,P) + Wj , ( Sl ) + ^(s 2 ,f) , h = !,-••,M (8) 


with W b (si) and V b (s 2 , f) are given by Equation (|5]), and Equation Q, respectively, 
and H° is the unbiased Hamiltonian for canonical molecular dynamics. For obtaining 
F(si,s 2 ), we require a strategy to obtain the unbiased probability distribution P(si,s 2 ). 

When no umbrella bias is present, under the quasi-stationary limit, the time depen¬ 
dent probability distribution from a WT-MTD simulation can be written as. 


P(R,0 = 


exp[-p{U(R) + V b ( S (R),0}] 


/dR exp [- j 6{U(R) +V b (s(R),f)}] 
and is related to the unbiased probability distribution P(R) by [17] 


(9) 


P(R) = P(R, f) exp 


P 


l/ b ( S (R),()-c(f)} 


( 10 ) 


where 


c(f) = - In 


f ds exp[—/3F(s)] 


( 11 ) 


J dsexp[—/3{F(s) + V b (s, f)}] 

Evaluation of c(t), however, requires a time-independent free energy F(s), which can 


be computed using the Tiwary-Parrinello time-independent free energy estimator [18] 


as. 


F(s) = — aV ,b (s, t) + — In / dsexp a.fiV h (s,t) 

P J l 


( 12 ) 
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where the time dependency of the first term cancels with that of the second. 

In WS-MTD, we construct the time independent probability distribution P^(si,s 2 ), 
by reweighting the MTD potential, for each umbrella h. On reaching the quasi-stationary 
limit, the statistical weight at which S 2 (R(f)) is sampled during MTD (due to MTD bias) 
is exp[—/3 {V^(s 2 (R),t) — C/,(f)}] (see Equation (10)). Thus, 

lt" dr ex p[£ W( s 2 (t),t) - c h ( r)}] S( Sl (r) - s\) 6{s 2 {t) - s' 2 ) 


Phi. s V s 2 ) “ 


(13) 


ft™* dr exp[/3 {VJ(s 2 (t),t) - c h ( t)}] 

In practice, for each h, we compute the above equation by discrete sum over the MTD 
trajectory of {s\(t),s 2 (t)}. The numerator is computed for bins (si,S 2 ) spanned within 
a chosen range, while the denominator is independent of the bin value. We compute 
integrals for a time series from f m i n to f max for which the quasi-stationary limit is 
applicable and a proper sampling of S 2 is obtained. In WT-MTD, quasi-stationary 
limit can be thought to be achieved when bias is growing slowly and uniformly in the 


domain of s 2 of our interest. The bias divergence law [28[ proves that C/,(f) oc ln(f) 
under this limit and thus it is preferred that the reweighting is carried out when this 
linear relationship is obeyed. |l8] 

It may be noted that P^(s\,s 2 ) is not reweighted for the umbrella bias w£(si). It 
is now straightforward to reweight P^(si,s 2 ) for the umbrella potential and combine 
this with WHAM such that slices of probability densities P^(si,s 2 ), h = 1, • • • ,M can 
be joined to obtain the unbiased distribution P(s\,s 2 ). WHAM involves minimizing 
the error in patching the M independent distributions {P^(si,s 2 )} by self-consistently 
solving the WHAM equations 


P(si,s 2 ) = 


Ltti n h exp \Pfh] exp[-j8Wfc(si)] 


(14) 


and 


ex p[—M = I ds 1 ds 2 exp[- J 6Wj , (si)]P(si,s 2 ) (15) 

where ft/, is the number of configurations sampled in the /z th window of the umbrella 
potential. The only difference with the usual WHAM equations is that biasing potential 
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along S 2 is set to zero. It is worth noting that Equation ( [14] > assumes that all the windows 
have nearly the same correlation time. |12|13] We have tested this assumption (following 
the procedure by Hub et al. 1291) for the severe case studied in Section |3.3| and found 
that this assumption is a valid one; see Supporting Information for details. 

Finally, the free energy surface F(si,S 2 ) is constructed using 


f(si,s 2 ) = -^lnP(si,s 2 ) (16) 

Thus, we have constructed a two-dimensional free energy surface by M independent 
1-dimensional WT-MTD with umbrella restraints. Although, here we have shown the 
WS-MTD equations for a two-dimensional case, it is straightforward to generalize this 
for higher dimensions. 

WS-MTD has the computational advantage that sampling can be parallelized over 
the M umbrella windows. However, the total computational time for WT-MTD in¬ 
creases with M and f max (Equation (fl3]>). The total simulation time required for each 
window, f max , depends on the time required to observe C/,(f) a In (f) behavior and the 
time required to sample CVs in MTD after reaching the quasi-stationary limit. We 
stress that the benefit of WS-MTD is mainly for the cases where a controlled sampling 
of a CV is required and cannot be achieved in standard MTD. Especially for the cases 
shown in Figure [lj, WS-MTD will be advantageous over WT-MTD. 


2.2 Efficient reweighting of near transition state regions in the MTD 
CV space 

To carry out reweighting, f max has to be large enough such that the bias potential 
is changing slowly and uniformly in the CV-space of our interest. If one needs to 
sample two minima separated by a large barrier within the MTD CV space (for a given 
umbrella), the first trajectory that crosses from the reactant basin to the product basin 
is insufficient for reweighting the regions near the transition state in the CV space, as 
the sampling near the transition state region is poor and has V h (s,t) » 0, thus quasi- 
stationary approximation is not applicable. On the other hand, for other parts of the 
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CV space where the system has visited several times, bias potential changes slowly 
and uniformly such that reweighting can be carried out. Ideally, we have to carry 
out long MTD simulation (for all relevant umbrellas) till the MTD trajectory recrosses 
the two minima multiple times so that bias potential grow slowly near the transition 
state region. However, this is often not practical, especially in ab initio simulations 
where the computational overhead for simulating recrossing trajectories is very high, 
and for many cases different CVs are required for the reverse reaction. Thus, when 
one is interested only in the forward process (and in computing the forward barrier), 
it is preferred that reweighting is limited to the reactant basin and near transition state 
regions in the CV space, but not the product basin. 

We use a simple and straightforward approach to reweight the transition state re¬ 
gions without simulating the recrossing trajectories. If a transition is observed for the 
first time from reactant well to the product well in the MTD CV space, we restart a 
WT-MTD simulation from an arbitrarily chosen point in the reactant well but using 
all the bias potential V b (s, t) accumulated in the previous WT-MTD simulation. When 
the reactant to the product transition is observed again, we repeat the same proce¬ 
dure. In this way, we increase the sampling near the transition state region, and the 
bias growth rate near the transition state region exponentially decreases towards zero, 
thereby reweighting can be carried out for the transition state region. Iterations are 
continued till a satisfactory convergence is achieved for the forward free energy barrier. 
This simple procedure has much less computational overhead and is used in the sim¬ 
ulations presented here. In practice, we always start the iterative procedure with the 
initial structure of the simulation, with velocities reassigned from Maxwell-Boltzmann 
distribution. 

2.3 Algorithm 

Here we briefly explain the algorithm for the WS-MTD approach. This technique re¬ 
quires no special implementation in a MD code which can carry out WT-MTD sim- 
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ulation and restrained dynamics simultaneously for different CVs. The reweighting 
procedure works as a post processing. 


1. For the chosen range of values of s a CVs, place M restraining umbrella potentials. 
For every umbrella, carry out a WT-MTD simulation sampling the coordinates. 
From these simulations, obtain the time series of s a (f),s^(f), V^(s,f) for some 
regular intervals of MTD time t. 


2 . 

3. 

4. 


Compute Ffr(s) using Equation (12) for h 


Compute C/ ; (f) using Equation ( [IT] ) for h 


1 ,••• ,M. 
1, • • • , M. 


Plot C/,(f), and based on that choose a time range, f m i n and t max , for which C/,(f) oc 
ln(f), and a proper sampling of has been accomplished. 


5. For the time range, construct the MTD-unbiased distributions P h(s a ,sp), h - 
1, • • • , M using Equation ( |13| ). It is crucial that the bin widths are chosen small 
enough to sample the fluctuations of every umbrella window. 


6. Using WFIAM based on Equation (14) and Equation (15), reweight the umbrella 
potential as well as combine the M distribution functions to get P(s a , s^). Note 
that P£( s a , Sp) can be the input for any standard WFIAM programs, but by setting 
the bias along the s« coordinates to zero. 


7. Using Equation (\6), construct the free energy surface F(s a , Sa 


2.4 Computational Details 

2.4.1 Two Dimensional Model System 


For testing the method, we considered a two-dimensional model system which has a 
broad basin and a narrow basin, whose potential is defined as 


U{x,y) = £u ; °exp 


Z —1 


Cli 


(x — X 


o\2 


+ (y-y°) : 


(17) 
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Parameters for the potential are given in Table [lj and the plot of U(x,y) is shown in 
Figure [2j The two minima are labeled as A (broad minimum) and B (narrow minimum), 
and the barriers for A—>B and B—s>A are 9.8 kcal mol -1 and 9.3 kcal mol -1 , respectively. 

Two different sets of simulations were carried out: (a) WS-MTD simulation with 
x coordinate chosen as the CV (si) for umbrella sampling and the y coordinate as 
the CV (S 2 ) for one-dimensional WT-MTD sampling; (b) a two-dimensional WT-MTD 
simulation with CVs x, and y coordinates (si, and S 2 , respectively). Umbrella potentials 
were placed along si from —0.05 to 0.44 Bohr at intervals of 0.01 Bohr. Restraining 
potential K/, for all the umbrella potentials were 9.9653 xlO 3 kcal mol 1 Bohr -2 . 

Mass of the particle was taken as 50 a.m.u. and constant temperature simulation was 
carried out using Anderson thermostat at 300 K. In WS-MTD and two-dimensional 
WT-MTD simulations, time steps for integrating the equations of motion were chosen 
as 0.12 and 0.17 fs, respectively. The initial Gaussian height (w 0 ) is 0.63 kcal mol -1 and 
the hill width parameter Ss is 0.005 Bohr. AT parameter for the WT-MTD simulations 
was taken as 1200 K. 

Before starting a WS-MTD, we carried out equilibration for a particular umbrella 
potential (for 2-5 ps), without adding MTD bias. Initial structure for an umbrella 
(away from the equilibrium) was chosen from the equilibrated structure of the nearest 
umbrella potential. This strategy has been also used for all the other systems studied 
here. 

2.4.2 1,3-Butadiene to Cyclobutene Reaction 

1,3-Butadiene can exist in cis and trans form and by an electrocyclic reaction it forms cy¬ 
clobutene. Free energy surface has a broad reactant basin and thus is an ideal problem 
to demonstrate the efficiency of the WS-MTD method using ab initio MD. 

In order to characterize the cis and trans isomers on the free energy landscape and 
the formation of cyclobutene, we have chosen the following CVs (see also Figure [3ja): 
a) distance C 1 -C 4 , d [Q — C 4 ]; b) the difference in the distances C 1 -C 2 and C 2 -C 3 , 
Ad[Ci — C 2 — C 3 ]. In WS-MTD simulations, d[Ci-C 4 ] was sampled using US and 
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Ad[C\ — C 2 — C 3 ] was sampled using WT-MTD. Two-dimensional WT-MTD simula¬ 
tions were also carried out to sample both CVs simultaneously. WS-MTD and WT- 
MTD simulations were performed within the framework of ab initio MD using plane- 


wave Kohn-Sham density functional theory (DFT) employing the CPMD program. [30] 
PBE exchange correlation functional [311 with ultrasoft pseudopotential [321 was used 


for these calculations. A cutoff of 30 Ry was used for the plane-wave expansion of 


wavefunctions. Constant temperature Car-Parrinello [33 [ MD at 300 K was carried out 


using Nose-Hoover chain thermostats. [34| A time step of 0.096 fs was used to inte¬ 
grate the equations of motion. A mass of 600 a.u. was assigned to the orbital degrees 

o 

of freedom. System was taken in a cubic supercell of side length 15 A. Extended La- 
grangian MTD approach was used here with a harmonic coupling constant of 1.0 a.u. 
to restrain the CVs and collective coordinates, and a mass of 50.0 a.m.u. was assigned 
to all the CVs. Langevin thermostat with a friction coefficient of 0.001 a.u. was used to 
maintain the CV temperature to 300 K. In WT-MTD simulations, Wq = 0.6 kcal mol -1 
and 5s =0.05 a.u. were taken. MTD bias was updated every 19 fs. Two indepen¬ 
dent two-dimensional WT-MTD simulations were performed using AT = 3000 K and 

00 o 

25000 K. In US, windows were placed from 1.5 A to 3.9 A at an interval of 0.05 A with 
Kh = 1.57 x 10 3 kcal mol ” 1 Ar 1 . In the case of WS-MTD, AT = 3000 K was used while 
all the other MTD parameters were the same as in the case of the normal WT-MTD. 


2.4.3 Controlled Sampling of Ligand Exchange in a Pd-Allyl Alcohol Complex in 
Aqueous Solution 

Here we are interested to compute the free energy barrier for the reaction WA1—»WA2 
(Figure [4]a) in aqueous solution. The purpose of this simulation is to achieve a con¬ 
trolled sampling, in particular, to avoid sampling the rapid ligand exchanges of trans 
Cl with solvent molecules to form WA3. [26 ] To simulate this process we have cho¬ 
sen four CVs: a) coordination number Pd to all the oxygen atoms of water molecules, 
C[Pd — O w ]; b) coordination number of Pd to Cl atom that is trans to the coordinated 
allyl alcohol, C[Pd — Cleans]; c ) coordination number of Pd to the two C atoms of allyl 
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alcohol to which it is coordinated, C[Pd — C]; d) coordination number of Pd to the Cl' 
in solution, C[Pd — Cl]. All initio MD simulations were performed with the same tech¬ 
nical details as mentioned previously. A periodic cubic box having the side length of 

° 

12.0 A containing the WA1 complex, one Cl - , and 53 water molecules were taken. The 


equilibrated structure of this system was taken from our previous study in Ref. [26]. 

We carried out WS-MTD simulations by sampling the CVs C[Pd — O w ] and C[Pd — Cleans] 
by US, and the CVs C[Pd — C] and C[Pd — Cl] by WT-MTD. In these simulations, the 
umbrella windows along the CV C[Pd — Cleans] were placed at 0.86 and 0.74, with K/, 

251 kcal mol -1 and 502 kcal mol -1 , respectively. The position of the umbrella window 
along the CV C[Pd — O w ] was fixed at 1.24. The K/, for the umbrellas along C[Pd — O w ] 
was 21.3 kcal mol -1 . The position and /q, for the umbrella restraints were based on the 
distribution of these CVs at equilibrium in the WA1 state. 

We carried out the extended Lagrangian MTD where a harmonic coupling constant 
of 2.0 a.u. was used to restrain CVs and collective coordinates, and a CV mass of 
50.0 a.m.u was assigned. Temperature of the CV was maintained to 300± 200 K by a 
direct velocity scaling. MTD biasing potential was updated every 29 fs and the MTD 
parameters vo o = 0.62 kcal mol -1 , and 5s = 0.05 were taken. AT = 3000 K was chosen 
for these simulations. 

As we are only interested in the forward barrier along the MTD CVs, we have 


used the strategy presented in Section 2.2 Moreover, to trace the minimum energy 


pathway on the reconstructed five-dimensional free energy surface, we used the string 
method. 


2.4.4 Free Energy for Water Coordination to the Active Site of NDM-1 by QM/MM 
Simulations 

To sample the reaction from Ell—>EI2 two CVs were chosen (Figure |5ja): a) distance 


between Zn \ and to the nearest water molecule (as identified in our earlier work [27]), 
d[Zrii — O w ]; b) coordination number of Zni to the carboxylate oxygen atoms Oi and 
O 2 , C[Zni — (Oi,02)]. The first coordinate was chosen as the US CV as we are expecting 
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a broad and deep free energy basin along this coordinate, while the other coordinate 

o 

was chosen for the WS-MTD sampling. Umbrella potentials were placed from 7.4 A to 
1.9 A at an interval of 0.1 A with k p = 4.5 x 10 2 kcal mol -1 A -2 . 

Hybrid QM/MM simulations were performed using the CPMD/GROMOS inter¬ 


face [361 as implemented in the CPMD package. The equilibrated structure of Ell was 
taken from our previous work |27j; see (Figure [5]). Two Zn ions, meropenem drug, 
side chains of Hisi 20 / Hisi 22 / Aspi 24 , Hisi 89 , Cys208/ and His 2 so and the water molecule 
near the active site were treated quantum mechanically (QM). Rest of the protein and 
the solvent molecules were treated by molecular mechanics (MM), and were taken in 
a periodic box of size 69.0x63.2x66.8 A 3 . The whole protein, including the QM part, 
and the solvent molecules were free to move during the MD simulations. The QM 
part was treated using the plane-wave DFT with Ultrasoft pseudopotentials and with 
a plane-wave cutoff of 30 Ry. A QM box size of 27.5x27.8x22.8 A 3 was chosen. When¬ 
ever a chemical bond has to be cleaved to define partition between QM and MM parts, 
boundary atoms are capped using hydrogen atoms. Capping hydrogen atoms were 
introduced between and C 7 atoms of Hisi 20 / Hisi 22 / Hisi 29 , and His 25 o and C a and 

atoms of Aspi 24 , and Cys208- 

Car-Parrinello MD was carried out for the QM part. Time step for the integration 
of the equations of motion was 0.145 fs. A mass of 700 a.u. was assigned to the orbital 


degrees of freedom and Nose-Hoover chain thermostats [34] were used to perform 
NVT ensemble simulations at 300 K. 

Extended Lagrangian variant of MTD was used where harmonic coupling constant 
was taken as 2.0 a.u. and the CV mass was set to 50.0 a.m.u. CV temperature was 
maintained to 300 K by coupling to Langevin thermostat with a frictional coefficient 
of 0.001 a.u. Gaussian potentials were updated every 29 fs and the MTD parameters 
Wo = 0.62 kcal mol -1 , Ss = 0.05 and AT = 7500 K were taken. We followed the strategy 
described in Section |Z2| for reweighting the CV space near the transition state along the 
MTD CV. 





-16- 


3 Results and Discussion 


First we benchmark the accuracy of WS-MTD method by sampling a two-dimensional 
model system where the free energy barriers are exactly known (Section |3.1[ >. To demon¬ 
strate the efficiency of WS-MTD simulation over the WT-MTD simulation in sampling 
broad and deep free energy wells, we model the cyclization reaction of 1,3-butadiene 
using ah initio MD (Section |3.2[ ). In Section 3.3 we demonstrate a controlled sampling 
of a high-dimensional surface using WS-MTD, which is otherwise difficult using nor¬ 
mal MTD. Subsequently, we present an example where WS-MTD is used together with 
QM/MM methods to sample the broad free energy surface of an enzymatic reaction, 
for which normal MTD simulations have failed. 


3.1 Two Dimensional Model System 

We carried out a two-dimensional WT-MTD simulation starting from the minimum A 
using the CVs Si(= x) and S2(= y). Several re-crossings were observed in the WT- 
MTD simulation, while a proper convergence in the barriers was observed by the third 
recrossing at 86.7 ps. The converged free energy barriers for A—>B and B—)>A are 9.8 
and 9.3 kcal mol -1 , respectively, and are identical to the exact barriers from the po¬ 
tential energy surface (Figure |2^). Free energy surface from this simulation is shown 
in Figure |2|a. We then carried out WS-MTD simulations, by sampling Si and S 2 using 
US and one-dimensional WT-MTD simulation, respectively. During these simulations, 
c(f) values for all the 50 umbrella windows were computed according to Equation ( [IT] ). 
For all the umbrellas, c(f)-f plots were nearly identical, and one of the plots is shown 
in Figure |2ja. c(f) oc ln(f) behavior was observed for t > 0.1 ps for all the umbrella 


windows. Thus, for reweighting the MTD simulations using Equation (13), we tested 
various values of f mm and f max greater than 0.1 ps, and the resulting free energy sur¬ 
faces were analyzed for convergence in free energy barriers. For f mm = 0.1 ps and 
fmax = 0.6 ps, the reconstructed free energy surface was quite noisy, and the free en¬ 
ergy barriers are having an error upto 0.4 kcal mol -1 (Table |2]). On the other hand, for 
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f m ax = 1-4 ps, free energy barriers for A—>B and B—were 9.8 and 9.2 kcal mol -1 , 
respectively, which is in good agreement with the reference free energy barriers; see 
also Figure [2j:. The topology of the relevant parts of the free energy surface is also 
well reproduced. For higher f max values, the free energy barriers converge to the exact 
values (Table |2|. The same convergence behavior was also observed for l m j n > 0.1 ps. It 
is also interesting to note that f m ^ = 0.0 and f max = 1.4 also gave reasonably accurate 
free energy barriers; see also the figure in the supporting information where error in 
the free energy barriers as a function of cumulative simulation time is shown for both 
WS-MTD and WT-MTD. Thus, performances of WS-MTD and WT-MTD are nearly 
the same, but the advantage of WS-MTD will become more obvious in the realistic 
examples discussed later. In cases with large f max , parts of the surfaces with higher free 
energy are better explored; see Figure [2ji,f,g. However, higher free energy regions are 
not often interesting and f max could be chosen based on convergence in the computed 
free energy barriers in more complex realistic systems. 

3.2 1,3-Butadiene to Cyclobutene Reaction 

Here we carried out a WS-MTD simulation starting with a 1,3-butadiene structure 
(Figure |3^). The converged reconstructed free energy surface from this simulation is 
shown in (Figure [3]:), where three minima CB1, CB2, CB3 could be observed corre¬ 
sponding to trans- buta-l,3-diene, ds-buta-1,3-d iene and cyclobutene, respectively. In 
WS-MTD, the converged barriers for going from CB1—»CB2 and CB2—»CB1 are 5.3 and 

3.3 kcal mol -1 , respectively. The free energy barriers for CB2—»CB3 and CB3—»CB2 are 
37.4 and 31.4 kcal mol -1 , respectively (Figure |3]:,d). Convergence in all the free en¬ 
ergy barriers is achieved with t m \ n = 0.0 ps and f max = 5.0 ps (for all the umbrella 
potentials); see Supporting Information for a more detailed analysis. 

To compare the performance of a two-dimensional WT-MTD sampling of both CVs, 
we carried out two independent WT-MTD simulations with two different AT param¬ 
eters. In the first simulation, AT was taken as 3000 K, which is the same as that was 
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used for WS-MTD. In this case, even after 1.2 ns of the simulation, no crossing for 
CB2—»CB3 was seen (see Supporting Information). The failure of WT-MTD is due to 
the presence of broad and deep reactant well and small AT parameter. Thus a second 
WT-MTD simulation was then carried out with AT = 25000 K, where we were success¬ 
ful in seeing several re-crossings from CB2—»CB3 within 1.2 ns. Although free energy 
barriers for CB1—^CB2 and CB2—^CBl have converged systematically (due to the small 
barriers separating them), the free energy barriers for CB2—>CB3 and CB3—»CB2 are 
not converged to the extent in WS-MTD simulations even after 1.2 ns (Figure |3]d). Free 
energy barriers computed from WS-MTD are well converged within 5.0 ps simulation, 
but not for WT-MTD even after 1.2 ns. Since we used 49 windows and f max = 5 ps 
per umbrella, the total computational time invested for a well converged result in WT- 
MTD is 245 ps. Also, note that all the 49 windows were running in parallel in our 
computations and the clock time to complete a converged WS-MTD was equivalent to 
the clock time required required for simulating 5.0 ps trajectory of a two-dimensional 
WT-MTD on same number of processors allocated per umbrella. 

We thus clearly show for a realistic system, that WS-MTD has better performance 
over the normal WT-MTD in sampling free energy surfaces with broad and deep basins. 


3.3 Controlled Sampling of Ligand Exchange in a Pd-Allyl Alcohol 
Complex in Aqueous Solution 


In our earlier work on the investigation of the Wacker oxidation of allyl alcohol [26] 
in aqueous solution we were interested in modeling WA1—>WA2. Flowever, we were 
unable to simulate this step using normal MTD due to rapid ligand exchange reac¬ 
tion of the trans- Cl (Cleans) with water molecules from solution, as a result of the 
strong frans-directing nature of the olefinic group (Figure [4^). Different wall potentials 
at different values of Pd-Cltrans coordination number were resulting in different free 
energy barriers and thus a reliable estimation could not be made. Thus, our earlier 


work [26 [ has used the equilibrium constant for WA2=WA1 from experiment and the 
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barrier for WA2^WA1 computed from MTD to estimate the WA1—>WA2 barrier as 
~19 kcal mol -1 . 

In order to compute the free energy barrier for WA1—^WA2, we have sampled the 
C[Pd — C], C[Pd — Cl], C[Pd — Cltrans] and C[Pd — O w ]. The controlled sampling of the 
C[Pd — Cltrans] and C[Pd — O w ] CVs was achieved by using US while the other CVs 
were sampled by a two-dimensional WT-MTD. The restraining values of US coordi¬ 
nates were determined form equilibrium simulation, and the restraining potential was 
computed from the standard deviation of the probability distribution of these CVs 
without any MTD bias. Along the C[Pd — O w ] coordinate, we do not want to sample 
other umbrella windows in order to prevent water coordination and thus the reaction 


to WA3. Based on our chemical intuition and from our previous experience, [26 [ we 
interpret that C[Pd — Cltrans] could vary from the equilibrium value along the reaction, 
because when external Cl - approaches axially and coordinate to Pd, an increase in the 
Pd—Cleans distance is expected. Thus we sampled another overlapping window for the 


C[Pd — Cltrans] CV at a lower value (Section 2.4.31. In this way, we carried out WS- 
MTD simulation to sample only the crucial part of the five-dimensional free energy 
landscape. Minimum energy pathway on this high-dimensional landscape was then 
obtained by the string method (Figure [4]:). [351 

Since reactant and product basins are present in the MTD CV space (see also Fig¬ 
ure |4]d) and that we are only interested in obtaining the forward barrier i.e. for the 
reaction WA1—»WA2, we used the iterative approach for reweighting as discussed in 
Section |2.2[ The convergence in free energy barriers with number of iterations is given 
in Table [3] A reasonable convergence was obtained for the third iteration, and free 
energy barrier for WA1—>WA2 is computed as 20 kcal mol -1 (Figure Up). This is in ex¬ 


cellent agreement with the estimated barrier of 19 kcal moP 1 in our earlier work. [26 [ 
By presenting this non-trivial example, we demonstrate that WS-MTD can be ap¬ 
plied for a controlled sampling of a high-dimensional free energy landscape of a com¬ 
plex reaction. Moreover, the iterative scheme for reweighting the near transition state 
regions is demonstrated for a realistic system. 
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3.4 Free Energy for Water Coordination to the Active Site of NDM-1 
by QM/MM Simulations 

In this section, we study a problem which is part of our ongoing research in elucidating 


the mechanistic details of antibiotic resistance in NDM-1. [271 In the intermediate stage 
of one of the proposed mechanisms for the hydrolysis of the acyl-enzyme complex 
within the active site requires a water molecule to enter the active site and coordinate to 
the Zni site, further leading to the dissociation of drug-carboxylate group coordination 
to Zip and subsequent protonation of /3-lactam N (see Figure [5jri) . |27| We attempted 
nearly 10 different MTD simulations with various collective variables to model the entry 
of a QM water to the active site and its coordination to Zn \, however, they all failed to 
simulate the reaction, and hinted a broad reactant basin. 

Thus, we addressed this problem using WS-MTD by US along d [Znj — O w ] and 
WT-MTD sampling along C[Zni — (0i,02)]. The US windows at d[ Zni — O w ] > 2.5 A 
was only simulated till the V h (s,t) ~ 0 limit is achieved in the relevant values of the 
C[Zni — (0i,0 2 )] coordinate. MTD biases are added for the windows with d[Zn\ — O w ] < 
2.5 A till a transition is observed along C[Zni — (0i,0 2 )] or till the total bias is equal to 
25 kcal mob 1 along the C[Znj — (0i,0 2 )] coordinate. If a transition is observed along 
the C[Zni — (0i,0 2 )] coordinate, we carried out the iterative approach for reweighting 
the transition state regions of the CV space as mentioned in Section 2.2| 

The reconstructed free energy surface thus obtained is given in Figure |5j:. The con¬ 
verged free energy barrier for this reaction is 23.6 kcal mol -1 . From the topology of the 
surface, it is now clear why the normal MTD simulations failed: the free energy basin is 
deep and flat along the d[ Zni — O w ] coordinate, and the minimum of the reactant basin 
is at d[ Zni — O w ] ~ 7 A. Free energy surface also indicates that water coordination to 
Zni and dissociation of Zni-carboxylate bond occur in tandem. 
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4 Conclusions 

Based on a combination of Tiwary-Parrinello and US reweighting techniques within 
WHAM, we propose the WS-MTD strategy for reconstructing complex free energy 
landscapes by sampling them by slices. Having the advantage of controlling the sam¬ 
pling along a few selected coordinates, WS-MTD aids in efficient sampling of flat, broad 
and unbound free energy wells with hidden orthogonal barriers. We also propose a 
simple strategy to reweight MTD trajectories for the regions near the transition state 
along the MTD CVs in order to obtain converged forward barriers without the need of 
sampling recrossing trajectories. The WS-MTD also has the advantage of paralleliza¬ 
tion over umbrellas. This feature enables sampling of high-dimensional landscapes in 
affordable clock time employing more computational resources in parallel. WS-MTD 
does not require any special implementation in MD codes that can simultaneously per¬ 
form restrained dynamics and WT-MTD. 

We carried out careful study on the accuracy of the method by sampling a two 
dimensional model potential where the free energy barriers are exactly known. Free 
energy barriers were found to converge systematically to the exact values with increas¬ 
ing f m ax- Several applications using the WS-MTD method were then presented. Ef¬ 
ficiency of WS-MTD in sampling broad and flat free energy surface is demonstrated 
by studying the cyclization of 1,3-butadiene using ab initio MD. WS-MTD yields con¬ 
verged free energy barriers with a much less computational cost (by a factor of 5 at least 
for this problem) than a normal MTD. By sampling the five-dimensional free energy 
surface for a ligand exchange reaction in a Pd-allyl alcohol complex solvated in water, 
we demonstrate the application of WS-MTD in controlled sampling of a high dimen¬ 
sional surface. The modeled ligand exchange reaction was earlier reported to fail using 
normal MTD. Also we showed the application of an iterative approach to reweight the 
MTD trajectories near the transition state regions on the CV space without the need of 
simulating multiple recrossing trajectories. Finally, we apply WS-MTD to sample an 
enzymatic reaction using QM/MM technique, where the coordination of a catalytic wa- 
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ter to the active site of NDM-1 enzyme is simulated. For this problem, WS-MTD could 
efficiently sample broad and deep free energy basins, which was otherwise unable to 
achieve in a normal MTD after several attempts. 

Since the method is practically applicable in ab initio MD and QM/MM MD simula¬ 
tions, we believe that this would be useful for sampling high-dimensional free energy 
landscape of complex chemical reactions. The technique would be beneficial for various 
problems in chemistry and biology such as A+B type reactions in weakly bound molec¬ 
ular complexes, reactions in Michaelis complexes, and substrate binding in enzymes to 
name a few, where one encounters broad and deep free energy basins. 
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Table 1: Parameters for the two-dimensional model potential 


i 

U° (kcal mol 1 ) 

d[ (Bohr 2 ) 

x° (Bohr) 

y°i (Bohr) 

1 

-17.885 

139.2985 

0.0000 

0.00000 

2 

-11.625 

41.7895 

0.3642 

0.00000 

3 

-11.625 

41.7895 

0.4916 

0.00000 

4 

-11.625 

41.7895 

0.6189 

0.00000 

5 

-17.885 

13.9298 

0.1821 

0.00000 
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Table 2: Free energy convergence in WS-MTD simulations for the double-well poten¬ 
tial. The converged free energy barriers from 86.7 ps of two-dimensional WT- 
MTD and the exact free energy values from the potential are also provided for 
comparison. 


Method 

train (ps) 

fmax (PS) 

A pi (kcal mol 1 ) 




A—^ B 

B A 

WS-MTD 

0.0 

0.6 

10.5 

9.5 



1.4 

10.0 

9.2 



2.8 

9.9 

9.3 


0.1 

0.6 

10.1 

8.8 



1.4 

9.8 

9.2 



2.8 

9.8 

9.3 


0.2 

1.4 

9.7 

9.1 



2.8 

9.8 

9.3 


1.5 

10.0 

9.8 

9.3 

WT-MTD 

- 

- 

9.8 

9.3 

Exact 

- 

- 

9.8 

9.3 
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Table 3: Convergence in the computed forward free energy barriers for WA1—»WA2 
with number of iterations. 


Iteration 

A pi (kcal mol x ) 

1 

19.4 

2 

19.9 

3 

20.1 




Free energy _ _ Free energy 
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Figure 1: Sketch of free energy surfaces showing three scenarios where traditional MTD 
simulations fail to sample efficiently. Green color shows the bias added by 
MTD. For cases (a) and (b), transition from reactant well (broad well at the 
right side) to the product well (narrow well at the left) would be failed to 
observe in a standard MTD. In the case of (c), computational overhead would 
be very high due to the broad reactant basin. 
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Figure 2: Two-dimensional double well potential: (a) Potential energy surface; (b) free 
energy surface from two-dimensional WT-MTD; Free energy surfaces from 
WS-MTD simulations for (t m \ n , t max ) values (in ps): (c) (0.1,1.4), (d) (0.1,2.8), 
(e) (0.2,1.4), (f) (0.2,2.8), (g) (1.4,11.4); (h) c(f) plot for one of the windows 
in the WS-MTD simulation where t axis is in log-scale. Free energy values 
are in kcal mol -1 relative to the free energy of the broad minimum (A). c(t ) 
is in kcal mol -1 . Contour values are drawn from 1.0 to 12.0 kcal mol 1 at 
1 kcal mol -1 intervals. 
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(a) 



c 2 c 3 


c 2 —C3 
C-( C 4 


CB1 CB2 CB3 



d [C-, -C 4 ] d [C-, -C 4 ] 



Figure 3: (a) Structure of trans- 1,3-butadiene (CB1), czs-l,3-butadiene (CB2), and cy¬ 
clobutene (CB3). Reconstructed free energy surface for 1,3-butadiene 
to cyclobutene reaction computed from (b) two-dimensional WT-MTD 
(fMTD=l-2 ns), and (c) WS-MTD (f m i n = 0.0, t max =15 ps). Free energy values 
are in kcal mol -1 relative to the free energy of the minimum (CB1). Contour 
values are shown from 1.0 to 45.0 kcal mol -1 for every 2 kcal mol -1 . CVs 

o 

are in A. (d) Various free energy barriers (color code: WS-MTD (black), WT- 
MTD (red)) on the free energy landscape computed for different f max values 
from WS-MTD simulation (bottom axis) together with the corresponding free 
energy barriers from two-dimensional WT-MTD with increasing MTD time, 
Alio (top axis). Here symbols ■, •, ♦, A represent free energy barriers for 
CB2—^CB3, CB3—s>CB2, CB1—^CB2, and CB2—>CB1 respectively (e) c(f) — t 
plot for one of the windows in the WS-MTD simulation, where the t axis is 
in log-scale. 
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5 1 

C [Pd-C] 


Figure 4: (a) Structure of Pd complex with allyl alcohol (WA1), PdCl^ (WA2) and Pd 
complex with water (WA3). Free energy barriers in kcal moP 1 are indicated 


over arrows, as computed from our previous study [26] (b)Snapshots from the 
simulation showing the complex of Pd with allyl alcohol (WA1), and PdCkj - 
(WA2); color code: Pd (tan). Cl (green), C (black), O(red), Fl(white) and sol¬ 
vent water molecules in blue stick representation, (c) Minimum energy path¬ 
way traced on the converged five-dimensional free energy landscape, (d) 
The non reweighted reconstructed free energy surface WA1—)>WA2 visual¬ 
ized along the MTD CVs for C[Pd — Cleans] = 0-74 and C[Pd — O w ] = 1.24. 
Free Energy values are in kcal mol -1 relative to the free energy of the mini¬ 
mum (WA1). Contour values are drawn from 1.0 to 20.0 kcal mol -1 for every 
1 kcal mol -1 . 
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Figure 5: (a) Structures of (Ell) and (EI2) with active site water molecule shown in blue 
color, (b) Snapshots of Ell, and EI2; color code: Zn (gray), C (black), O (red), 
N (blue), S (yellow), and H (white). Active site residues and the attacking 
water molecule are in CPK representation, while the drug molecule is shown 
in stick representation. (c)Reconstructed free energy computed from WS- 
MTD. Free energy values are in kcal mol 1 relative to the free energy of the 
minimum (Ell). Contour values are drawn from 1 to 26 kcal mol -1 for every 
2 kcal mol -1 . d[ Zni-O w ] is in A. 
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Abstract for the Table of Content Figure: 

Here we present a technique to sample a high-dimensional free energy landscape 
as slices by combining umbrella sampling and metadynamics simulation sampling or¬ 
thogonal coordinates simultaneously The full free energy surface is then reconstructed 
by combining these slices using the Weighted Histogram Analysis Method. 






